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We present a semi-analytical free-energy model aimed at characterizing the thermodynamic prop- 
erties of dense fluid helium, from the low-density atomic phase to the high-density fully ionized 
regime. The model is based on a free-energy minimization method and includes various different 
contributions representative of the correlations between atomic and ionic species and electrons. This 
model allows the computation of the thermodynamic properties of dense helium over an extended 
range of density and temperature and leads to the computation of the phase diagram of dense fluid 
helium, with its various temperature and pressure ionization contours. One of the predictions of 
the model is that pressure ionization occurs abruptly at p 10 g cm -3 , i.e. P J> 20 Mbar, from 
atomic helium He to fully ionized helium He 2+ , or at least to a strongly ionized state, without He + 
stage, except at high enough temperature for temperature ionization to become dominant. These 
predictions and this phase diagram provide a guide for future dynamical experiments or numeri- 
cal first-principle calculations aimed at studying the properties of helium at very high density, in 
particular its metallization. Indeed, the characterization of the helium phase diagram bears im- 
portant consequences for the thermodynamic, magnetic and transport properties of cool and dense 
astrophysical objects, among which the solar and the numerous recently discovered extrasolar giant 
planets. 

PACS numbers: 52.25.Jm, 05.70.Ce, 52.25.Kn 



I. INTRODUCTION 



Within the past decade, over a hundred brown dwarfs, astrophysical bodies not dense enough to sustain hydrogen 
fusion in their core, and extrasolar giant planets, i.e. jovian planets orbiting stars outside the solar system, have 
been discovered. These objects are composed essentially of hydrogen and helium. Given their large gravity and 
relatively low temperature, within astrophysical standards, the hydrogen and heluim fluid is under an atomic or 
molecular form in the outermost part of the body and under the form of a fully ionized electron-ion plasma in the 
innermost regions. Such an internal structure is common to many so-called compact objects, from our own jovian 
planets to the external layers of white dwarfs or neutron stars. The characterization of the structure and cooling 
properties of these compact objects thus requires the knowledge of the thermodynamic properties of dense hydrogen 
and helium fluids, and more importantly a realistic description of the partial, pressure ionization regime. Given 
the large variations of thermodynamic conditions characteristic of the structure and evolution of such astrophysical 
bodies, these thermodynamic properties, characterized by the equation of state (EOS), must be calculated over several 
orders of magnitudes in density and temperature. As discussed below, the necessity to calculate the thermodynamic 
properties over such a large range of conditions precludes the use of heavy computer simulations and thus necessitates 
the derivation of EOS models which allow extensive calculations within a reasonable amount of computer time, 
unfortunately at the price of a more approximate, or say phenomenological description of the properties of matter at 
high-density. 

Interestingly enough, these EOS of dense matter under astrophysical conditions can now be probed on Earth by 
shock wave experiments. Future large laser experiments, like e.g. the NIF project at Livermore or the LMJ project 
in France, will reach conditions characteristic of the deep interior of the aforementioned astrophysical bodies. So 
not only the calculation of dense matter EOS is of interest for astrophysical applications, but it is necessary for the 
confrontation of theory with existing and future high-pressure experiments, yielding eventually a correct knowledge 
of the properties of matter under extreme conditions. Hydrogen, the most common element in the universe, has been 
studied extensively, both on the experimental and theoretical fronts, and the EOS of dense hydrogen becomes more 
and more constrained, although the very regime of pressure ionization still remains ill-determined. The same cannot 
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be said for helium. Although some experiments exist in the regime of neutral helium at high-density, as detailed below, 
the regime of helium pressure ionization, from He to He + and He 2+ remains for now unexplored, and no attempt has 
been made to give a detailed theoretical description of these domains. It is the very purpose of the present paper 
to derive an EOS for dense, partially ionized helium, covering the gap between the previous study of dense neutral 
helium jl| and the fully ionized regime |2| As mentioned above, not only the calculation of such a dense helium 
EOS is necessary for a description of the thermodynamic properties of astrophysical compact objects, in particular 
the recently discovered gaseous exoplanets, but it provides a useful guide for future high-pressure shock-wave or laser 
experiments. 

The paper is organised as follows. In Sec. ^ we briefly comment on the general formalism underlying the present 
calculations. The various contributions entering our general model free-energy are presented in detail in Sec. II I II 
The results, and the limitations of the model are presented in Sec. IIVI Special attention is devoted to the impact of 
various approximations in the free-energy calculation on the final results. Section is devoted to the conclusion. 



Equation-of-state calculations can be divided into two generic categories. The "physical approach" is formally 
exact as it involves only fundamental particles, electrons and nuclei, interacting through the Coulomb potential. The 
partition function is calculated using the eigenvalues corresponding to this TV-body system. In practice, however, 
the exact solution cannot be calculated, in particular when bound-states form, and either perturbative expansions 
or approximate numerical schemes must be used. The validity of the expansions is limited to high temperatures 
and/or low-densities, i.e. apply to weakly or moderately coupled plasmas. The regime of pressure ionization thus 
cannot be described by such expansion schemes. Numerical technics, such as density functional theory, molecular 
dynamics or path-integral Monte Carlo simulations, do extend to the strongly correlated regime but the description 
of the pressure ionization regime then becomes a formidable task, and involves also physical approximations in the 
calculations of either the electron functional or the nodal functions, not mentioning the finite size effects due to 
the limited number of particles in the simulation. In practice, these simulations do not allow the calculation of 
thermodynamic quantities over a large range of temperatures and densities, as needed for practical applications, 
as mentioned earlier. For this reason, a more phenomenological approach has been developped which combines a 
simplified description of the properties of dense matter and a semi-analytical derivation, allowing the calculations of 
extended thermodynamic tables with moderate computer time investment. This is the so-called "chemical picture" . 
In this approach, the basic particles are no longer only electrons and nuclei but also bound species (atoms, molecules, 
ions), which are characterized by their interparticle potentials. That means that the particles remain distinguishable 
(in a classical sense) in the plasma, with their own identities and interaction properties. The problem thus reduces to 
the free-energy minimization of a multi-component system, taking into account chemical and ionization equilibrium 
between the various species. Although certainly of doubtful validity in the regime of pressure ionization, where the 
concepts of pair potential and bound-states become meaningless, this approach has been shown to yield reasonably 
accurate descriptions of hydrogen at high density Moreover, as mentioned above, this approach presents the 
advantage of being semi-analytical and thus has a precious practical interest for EOS calculations. Last but not 
least, the chemical approach offers the noticeable advantage of clearly identifying the terms and the approximations 
aimed at describing various physical effects. Such terms can be added or removed with limited effort, allowing a rapid 
identification of the dominant contributions responsible for the thermodynamic properties of matter under complex 
conditions. Therefore, despite its shortcomings, the chemical approach should be seen as a useful alternative to the 
"exact" physical approach. 



The chemical approach is based on the minimization of the free-energy F({Ni},T,V) corresponding to a system 
containing {iVj} different species inside a volume V at temperature T. This minimization SF = -^-8Ni = must 
satisfy the electroneutrality condition and the stcechiometric conditions corresponding in our case to the following set 
of chemical equations: 



II. 



GENERAL CONSIDERATIONS 



A. Chemical picture of a dense plasma 



B. General free-energy model 




(1) 
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The canonical partition function of the system Z is assumed to be factorizable into different contributions, so that 
the free-energy F = —kT\nZ can be split into the sum of translational, configurational and internal contributions [(| 
■ Adding up the correction arising from the quantum behaviour of the heavy particles, one gets: 

F({Ni},T,V) = F id ({Ni},T,V) + FcorfdNi},^) 
+ F int ({N z }, T, V) + F qra ({Ni}, T, V). 

(2) 

The conditions of validity of such a separability are: 

• the discretization of the eigenvalues corresponding to the translation degrees of freedom and to the center-of-mass 
positions are negligible. This is the quasi-classic approximation; 

• there is no coupling between the translation degrees of freedom and the center-of-mass positions; 

• the internal energy levels remain essentially unperturbed by the interactions with surrounding particles. 

If the two first conditions are satisfied in the present context, the last one certainly becomes invalid in the pressure 
ionization regime. We expect this regime, however, to cover a limited range of density, as pressure ionization generally 
occurs rather abruptly. Eventually, only comparison with experimental data can give a quantitative estimate of the 
discrepancy due to this underlying factorization condition. The various contributions to F are described in the next 
section. 

III. FREE-ENERGY MODEL 

We first present the models used to calculate the contributions to the total free-energy arising from each different 
species, He, He + , He 2+ . Then, we describe the modelization of the coupling between these various species. 

A. Model for atomic helium He 

1. The kinetic free-energy _F|d 
The ideal part of the free-energy, corresponding to the kinetic part of the Hamiltonian, is given by 

V (2-KMk B T^ 3/2N 



F id (N, T, V) = -Nk B T 



1 + ln 



(3) 



N \ h 2 

where N is the number of helium atoms of mass M inside the volume V at temperature T. 

2. The configurational free-energy F co ^{ 

The configurational free-energy F con { , arising from the interactions between helium atoms, is calculated within the 
Weeks- Chandler- Andersen [2j [10] (WCA) perturbation theory. The interaction potential <&(r) is split into a reference 
potential < i ) r of(r) and a perturbative part ^pcrtM- Truncating the perturbative expansion of the free-energy after the 
first order, the so-called high-temperature approximation (HTA), yields: 

N 2 f 

F conf = F rcf (T, Vj N) + — J dr $ port (r)g rcf (r, V, N) . (4) 

The problem is thus reduced to the potential separation and to the calculation of F TC [(T,V, N) and g Te f (r, V, N). 
Concerning the first point, we use a modification of the procedure of Kang et al. pl|. namely: 

rcfl ' 10 if r > A, 

$ tM - J <KA) + f|,=A(r-A) if r<A, [b) 
pcrtl ' ~ { $(r) if r > A, 
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FIG. 1: Hard-sphere radii ctwca for atomic helium (solid line), for T = 10 3 K, T = 10 3 ' 5 K and T = 10 4 K from top to 
bottom, and break-point A (dotted line) as a function of the density. 



where A = (a^, 3 + r*~ 3 )~ 1 / 3 ; af cc = {s/2/ '(N /V)) 1 ^ 3 and r* corresponds to the minimum of the potential $(r). This 
choice for the density-dependent break-point A has the advantage to give a continuously differentiable A. Concerning 
the second point, we approximate the repulsive reference potential by a hard-sphere potential. The hard-sphere radius 
a is calculated from the Barker-Henderson criterion: 



0"BH = 



/ dr (l-e-P*™*) = / dr (1 - e^*"*), (6) 
Jo Jo 

with the Verlet & Weiss correction to include a density-dependence: 

CT = - H ( 1 + ^)' (?) 

where S is a function of the temperature and ct\/2(Jq is a function of T and a. This non-linear equation is solved by 
direct iteration, using obh as an initial guess for a = a(T,n). An example of the evolution of a with density and 
temperature is presented in Fig. ^ The free-energy and the radial distribution functions for the hard-sphere reference 
system F rc f = Fhs, Srcf = 3hs, are obtained analytically [l4| . 

To describe the interaction between two helium atoms, we choose the Aziz & Slaman [l5j potential for r > 1.8 A, 
and the Ceperley & Partridge [la] one for r < 1.8 A. Following Aparicio & Chabrier this two-body potential is 
modified by a density-dependent function to mimic the softening due to A-body effects at high density: 



*w=((i-c) + ^y^o 



(r). (8) 
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FIG. 2: Interaction potential between two helium atoms, without iV-body correction (solid line) and with the iV-body softening 
correction at p = 10 g/cm 3 (dotted line). 

The two parameters C and D are optimized to reproduce the experimental measures of adiabatic sound velocity 
[r7|. A x 2 minimization yields (C,D) = (0.44,0.8 cm 3 /g). This potential is illustrated in Fig. [21 whereas Fig. |3] 
compares the measured sound velocity and the one calculated with our potential. Fig. ^ compares the Hugoniot 
curves obtained with the present atomic helium free-energy model and interatomic potential with presently available 
shock- wave experiments [l8|. These comparisons assess the validity of the present model at least up to the limit of 
the data, i.e. P~ 1 Mbar. 

3. The internal free-energy fi nt 

The divergence of the internal partition function, gi exp(— Ei / ksT) , of an isolated atom is a well-known problem 
in statistical physics. It emphasises the necessity to take into account the interactions between atoms in the calculation 
of the internal partition function, Z; nt = exp(— (3F[ nt ). For a density n, each atom has a typical available volume n^ 1 / 3 
so that, as density increases, the levels associated with the highest eigenvalues will move into the continuum. When 
the density is high enough to disturb even the ground-state, the electrons can no longer remain bound to the nuclei: 
this is the pressure ionization phenomenon. We have included the effect of the interactions of surrounding particles on 
the internal partition function of helium within the so-called occupation probability formalism (OPF). The OPF 
ensures the statistical-mechanical consistency between the configurational free-energy characterizing the interactions 
between atoms, -Fconf, and the internal free-energy contribution, Fiat. The OPF has been extensively presented in 
various papers (see e.g. pofl. and is only briefly outlined for completeness. We consider a system of interacting 
particles, of free-energy F = F^ — fc^Tln Z lnt + /, where / is the non-ideal term. Within the OPF, the total 
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P (g/cm ) 

FIG. 3: Comparison between the experimental measures of the adiabatic sound velocity [l7j (solid line) and the present 
calculations with (C,D) = (0.44,0.8 cm 3 /g) (crosses) and (C,D) = (0,0 cm 3 /g) (dotted line). 



free-energy can be rewritten under the form: 



F = F id -k B TlnZ + f-J2 N a 



df 

dN a 



(9) 



with 



Z = y^ a ff a e" 



/3E a 



and ui a — exp — j3 



df 

dN a 



(10) 



The term u> a can be seen as the probability that the eigenstate a of the atom still exists in the midst of the surrounding 
particles. These factors u> a are calculated consistently from the configurational term /, and the term *}Z a N a Jf^- 
ensures the statistical-mechanical consistency (see [21 )• The OPF has several noticeable advantages, among which: 

• Lo a decreases monotoneously and continuously with increasing density, ensuring the convergence of Z- mt and the 
derivability of F int ; 

• no ill-controlled energy shifts of the levels are introduced , as required from the condition of factorizability of the 
partition function (Eq. |2J). Experiments at low-density |2l| and calculations [2^| do not show such energy 
shifts; 

• the probabilistic interpretation of u) a enables us to combine several occupation probabilities arising from statis- 
tically independent interactions. We will come back to this point in Sec. IIIIDI 
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FIG. 4: Comparison between the experimental single and double-shock Hugoniot curves [Tsl ] and the present calculations with 
(C,D) = (0.44,0.8 cm 3 /g) (solid line) and (C,D) = (0,0 cm 3 /g) (dotted line). 

The exact solution, in principle, requires the knowledge of all the interaction potentials between an atom in state a 
and an other one in state a' . In the absence of such information, we have adopted the simplest approach which consists 
to characterize excited state interactions by hard-sphere excluded volumes in the phase space. The hard-sphere radii 
are calculated with the scaling law derived by Aparicio & Chabrier (Eq. (14a) and (14b)). Within the first order 
in the expansion of the non- ideal part / of the free-energy (Eq. |2J, the ui^ s for the excited states are thus given by: 

This nonlinear equation is solved iteratively by using results obtained within the low excitation approximation (LEA) 
and low density approximation (LDA), w hs,lea+lda _ e xp(— 7tA^(ct q + <Ji) 3 /6V), as initial guess. 

4- The quantum correction of the free-energy F qm 

We have taken into account the correction to the free-energy arising from quantum effects due to the finite size of 
the atoms by keeping the first order of the Wigner-Kirkwood fi 2 expansion of Tr [e - ^] |24| [2{| 

F Qm = N 2 / dr V 2 $(r)o(r). (12) 

q 24fc s TV r M H c J 

^(r) corresponds to the potential explicited in Sec. IIII A 21 and g(r) = y{r)e~ '' 3 *( r ) is approximated by y-Hs{f)e~^ (j '\ 
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B. Model for the partially ionized plasma {He + , e } 

Because of the presence of bound-states, the treatment of He + presents the same difficulties as for He. We adopt the 
same formalism, namely the WCA perturbation expansion, to calculate the He + configurational free-energy (with the 
hard-sphere model as the reference system) and the OPF to treat the internal partition function. For the long-range 
interaction potential between He + ions, we take a Yukawa potential, e~ ksT /r, where the density- and temperature- 
dependent screening wave vector is given by |26j| : 

k s (n, T) = ^=kT F [6 x ' 2 F_ 1/2 ^/k B T)] 1 ' 2 , (13) 

where fcxF = (4:Tn e e 2 /nh 2 ) 1 ' 2 (3Tr 2 n e ) 1 ' 6 is the Thomas-Fermi screening wave vector, n e is the total free electron 
density, 9 = T/Tp is the electronic degeneracy parameter (Tp is the electron Fermi temperature), F n is the Fermi 
integral of order n, and /i/kgT is the electron chemical potential defined by Fi/^ip/kgT) = 2#~ 3 / 2 /3. 
For the treatment of the internal free-energy, we need a scaling law to associate a hard-sphere radius to the excited 
states of He + . Since He + is hydrogen- like, and the energy levels are degenerate toward the orbital quantum number 
I, we write this scaling law as: 

<r n = n 2 <7 1 , (14) 

where u\ is the WCA hard-sphere radius associated to the ground-state, and n is the main quantum number. 
The calculations then proceed exactly as in Sec. IIII Al 

C. Model for the fully ionized plasma {He 2+ , e } 

The free-energy of a fully ionized electron-ion plasma (FIP) has been calculated by Chabrier & Pothekin 0] and 
Potekhin & Chabrier These authors provide analytical parametrizations for the various thermodynamic quantities. 
We refer the reader to these papers for a description of the fully ionized plasma model. 

D. Interactions between different species 

Besides all the aforedescribed contributions to the free-energy, arising from interactions between species of same 
nature, we must also include contributions arising from the interactions between species of different nature. 

1. Hard-sphere interactions between atoms and ions 

The first order interaction between the atomic and ionic species He, He + and He 2+ is the hard-sphere excluded 
volume interaction, FHs({^He,a, N Re + a , iV He 2+}, {cthc,q, c Ho +,q, o" Ho 2 +}, V, T), with a radius ohc 2 + = f° r the He 2+ 
ions, calculated consistently from the hard-sphere free-energy of a multicomponent interacting system [l3j |. It can 
be shown easily that the contribution arising from the o"He 2 + = component is equivalent to renormalizing the 
ideal (kinetic) term for this species with a volume V = (1 — T])V, where i] = J2ie{ns} 7TUi ' T i corresponds to the 
total packing fraction p|. This term thus takes into account the He-He, He + -He + , He-He + , Hc-Hc 2+ and He + -He 2+ 
interactions. Note that, contrarily to previous approaches, we do not consider excluded volume interactions between 
bound species and free electrons. Indeed, such an approach does not seem to be justified, for the quantum exclusion 
principle applies only to electrons in the same state. The entire volume of the system is thus available to the majority 
of the free electrons, even in the presence of bound species, as far as the free electrons are in a quantum state different 
from those corresponding to the bound-states. In any event, we have checked that the introduction of an excluded 
volume for the electrons does not modify significantly the final results. 

2. Induced interactions between atoms and ions 

The presence of charges in the neighbourhood of species with bound-states has two consequences. The first one is 
the induced polarization due to the electronic cloud, which translates into a related contribution to the free-energy. 
The second one is the induced Stark effect on the bound-states, due to the ambient electric field which modifies 
the one associated to the atom nucleus. These two effects have been taken into account in our model as described below. 
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a. Polarization effects 

The polarization contribution to the free-energy arising from the interactions between the charges and the neutral 
atoms He has been handled as in previous TV-body approaches [22j 

F pol = —?-N Hc J2 N i B He,i- (15) 

i=Hc+,He 2 +,e- 

The second virial coefficients Bu c ,i are given by: 

poo 

B ne>i = 2n drr 2 (l-e-P*U), (16) 



where 

Z, : e 2 a,; 



*pol« = - 



1 + k s r 



'Hc-i 



2 

exp (-2k s r) (17) 



is the polarization potential between He and the species i. The two free parameters OHe-i and cti are the hard core 
radius and the polarizability. For the He-He 2+ and He-e~ interactions, the hard core radius is chosen to be the He 
atom ground-state radius, c^_ Hc , and the polarizability (which has the dimension of a volume) is equal to (er^_ He ) 3 . 
For the He-He + interaction, the hard core radius is c£p_ H + = (opfc-Hc + °h +-H +)/^ ano - the polarizability is equal 
to«-H C+ ) 3 - ' ' ' 

b. Electric microfield effects 

Stark effect on the bound-states, arising from the electric microfield E due to the surrounding charges, is also 
treated within the framework of the OPF. The occupation probability associated with the Stark interaction on the 
internal states of He and He + is given by 0] : 

ocrit 

= / " dpPtf), (18) 



where f3 = (ineoa 2 /Ze)E is the dimensionless electric field (Ze is the ion charge and a — (47rn/3) -1 / 3 is the mean 
interparticle distance), V((3) is the probability that the central ionic center experiences a field between j3 and /3 + d(3, 
and /S™* is a critical field associated to each bound-state a. Potekhin et al. 28} have calculated the microfield 
distribution of an atom (neutral ionic center) or an ion (charged ionic center) immerged in a surrounding ionized 
plasma. These calculations take into account the interactions in the plasma (r = (Ze) 2 /akT ^ 0), and recover the 
Holtzmark limit in the case of a non-interacting, perfect gas (r = 0). These authors provide analytical formulae 

for Q(P,Y) = fgdtP(t,T) in the case of a neutral or a charged central ionic center. Note that Q(j3,T) and thus 
the probability cj^ e not only depend on the temperature, as in the Holtzmark limit, but depend also on the density, 
through the parameter V. The critical fields are given by Hummer et Mihalas in the case of a hydrogen-like system. 
We have directly applied their prescription to He + , and used the similarity between a He atom and a hydrogen-like 
system, with a central charge equal to 7/4 for the ground-state and 1 for the lsnl-type levels Q, to calculate the 
critical fields corresponding to atomic helium He. 

3. Long range interaction between He + and He 2+ 

The remaining coupling contribution between the various species stems from the long range Coulomb interaction 
between helium ions He + , He 2+ and electrons. Short distance interactions due to the internal levels of He + have 
been considered in the previous sections. The treatment of the long range Coulomb interaction between the two ionic 
species will certainly have some impact in the pressure ionization regime where He + and He 2+ coexist, but will not 
modify the rest of the phase diagram. This contribution, however, is difficult to evaluate accurately. Considering the 
He + -He 2+ interaction as a pure Coulomb contribution, thus representing the He + -He 2+ fluid as an interacting two- 
component Z\ = 1, Zi = 2 point-charge plasma is not satisfactory, for it precludes a correct treatment of the internal 
levels of He + , which has been included in our formalism (see previous section). In this context, and in the absence 
of an accurate formulation, we estimate the contribution to the total free-energy arising from the He + , He 2+ , e~ long 
range interaction in the framework of the ion-sphere model |29j , thus considering only the electrostatic contribution 
to the free-energy. In this very simplified model, the interaction between He + and He 2+ gives a contribution equal to 
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ZiZ-ie 1 1 a per pair, with Z\ — 1 and Z^ = 2, whereas the contribution due to the interaction between the central ion 
He + and the uniformly charged sphere —Z-ye gives a contribution — 3/2(Z 1 e) 2 /a per He + . The He 2+ -e~ and e~-e~ 
contributions are already included in the FIP model mentioned in Sec. IIII CI The contribution thus reads: 

F(N N \ ^Hc+Wh^ e 2 3 (Z ie ) 2 

^(A He +, 7V He 2+J = o Z 1 Z 2 oJ V He+ ■ U 9 ) 

The very crude treatment of this interaction between He + , He 2+ and electrons is certainly a major shortcomings 
of the present model and Eq. 1191 gives at best an order of magnitude of the contribution of this interaction to the 
free-energy. As mentioned above, there is no satisfactory description of ions with bound-states, He + in the present 
context, immersed in a surrounding dense plasma. Indeed, it is difficult to capture the drastically different nature of 
the short-range and long-range interactions of such species with surrounding charged particles. This is undoubtedly 
a limitation of the chemical picture, and of the related distinction between different entities. In reality, the concept 
of identifying He + or He 2+ particles, based on a concept of potential or pseudopotential, becomes meaningless at 
high density. Only at high temperature, when kinetic contributions dominate, is the approach conceptually correct. 
Therefore, although Hc + or He 2+ are distinguishable in our model free-energy, we do not pretend to give an accurate 
description of the second stage of helium pressure ionization, from He + to He 2+ . As detailed in the next section, 
however, we have checked that the present, crude description of the He + -He 2+ interactions does not alter the final 
phase diagram. The reason is that, at least in the present model, helium pressure ionization proceeds directly from 
atomic helium He to fully ionized helium He 2+ , or at least to a strongly ionized stage. It will certainly be interesting 
to compare these results with experiments and with results obtained with first-principle calculations, although these 
latter will certainly have to face their own difficulties in this complex regime. 



4- Summary 

Summarising out the various contributions described in the previous sections, and following Eq. 0] and [§1 the full 
model free-energy reads: 



(3F 



(V,T,{N i }) = - £ 



=Hc,Ho+ 



A/to, 



, / V {2nMk B T 
1 + In — 



Ni 



h 2 



3/2 N 



=Ho,Hc+ 



a w IiS w ME e"' 3£; " 



PF KS ({iV H e,a , ^He+ ,c - },Wiie = 0},V,T) 



+ 



y 

i={Hc,a},{Hc+,a} * 0t 
i j— He, He" 1 " 



Ni pdF m ({N l },{a iio , a ,c7 lic+ ^},V,T) (3F po \ (He, He + , He ,N e -,V,T) 



ON 



N u 



N t , 



24(fc B T) 2 M Ho V ^ N tot N tot 

i, j— He, He" 1 " 



J dr V 2 $(r)e- /3 * M Wyg s (r 



+ 



Nn e +N Be ?+ 



Z1Z2 



3N Re+ (Z ie y , /3^ iF (y,T,7V Hc2+ ,iV e -) 



2A^ ot a 2 Not a N tot 



(20) 



where N tot = N He +N Ue + + N He 2+ . Note that 

UJ ia' UJ ia include the occupation probabilities calculated from interactions 
with neutral surrounding particles (hard-sphere interaction, Eq. 11 If) and with charged surrounding particles (microfield 
interaction, Eq. 118)) . 

The equilibrium populations are derived from the minimization of the frcc-cnergy F(V, T, {N}) with respect to two 
independent variables, given the conditions of mass conservation, A^ He 2+ = A^ot — A^ Ho + — A^je, and electroneutrality, 
A^ e - = A Hc+ +2A Hc 2 + : 



U N tot 



cW Hc 







d 



0F 
N to% 



(21) 



T,V,JV H o 



dN He+ 

T,V,N Ho+ 

Convergence of this two-dimensional minimization is achieved when the change in the populations from one iteration 
to the next one is less than one part in 3 10~ 7 . The various thermodynamic quantities are then calculated from 
appropriate derivations of the free-energy. 
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FIG. 5: Comparison between the populations obtained with our model (lines) and those corresponding to the Saha equations 
(symbols). The solid (resp. dotted) line corresponds to the He (resp. He + ) fraction. The temperature is T = 10 4 ' 7 K. 

IV. RESULTS 

As mentioned previously, our free-energy model, with the He-He potential calibrated on sound velocity measure- 
ments reproduces the available Hugoniot experiments ^| (see Fig. 0J. We have also checked that we recover 
the results of the Saha equations in the low-density limit and, by construction, the fully ionized plasma model at high 
density. An example is shown in Fig. [5] for T = 10 4 ' 7 K. The vanishing fraction of bound species populations for 
p ^ 10 g cm~ 3 illustrates the onset of pressure ionization. 

We have also checked that we recover the results of Aparicio & Chabrier for pure atomic helium in the low-density, 
low-temperature regime until pressure ionization sets in (see Fig. [5]for T=10 3 ' 5 K). 

A. Limitations of the model 

As mentioned earlier, our free-energy minimization method is rooted in the chemical approach. It is based on a 
heuristic treatment of the dominant physical effects responsible for the thermodynamic properties of dense atomic 
or ionized helium. Although it certainly retains some degree of reality, this model cannot pretend giving an exact 
description of these properties, and the results should depend to some extent on the main approximations used in the 
model. We examine this issue in the present section. 
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FIG. 6: Comparison between the present model (solid lines) and the results of Aparicio & Chabrier 0] (dotted lines with 
symbols), which do not include pressure ionization, for the pressure, the massic entropy, the massic internal energy and the 
specific heat at T=10 3 ' 5 K. 

1. Lower boundary for o\ 

At very high density, the WCA radii tend eventually to zero, as shown in Fig. 2] This favors the He and He + 
species and thus prevent pressure ionization to occur, a well identified artifact of the chemical picture |^| |3(| . In 
order to prevent such an unphysical behaviour, we define arbitrarily a lower limit for <7He and cr Hc +. Fig. (for the 
10 4 7 K isotherm) illustrates the effect of this approximation for o\ > 0.8 A and <j\ > 0.5 A. 

Not surprisingly, the choice of a lower limit for <j\ affects appreciably the populations in the very regime of pressure 
ionization. However the effect is almost inconsequential on the thermodynamic quantities, the very purpose of the 
present calculations. This stems from the fact that the bound-species do not contribute to the free-energy when they 
are associated with a very small radius. The final model calculations were made with o\ — 0.8 A. 

2. Polarizability of He-He + 

We have also tested the influence of the polarizability a,i which appears in the He-He + potential, and which has 
been taken equal to the volume o"hc-Hc+ ■ Calculations conducted with a value of a; reduced or increased by a factor 
10 left the results nearly unaffected. This can be easily understood as in the domain where non-ideal effects play a 
role, He and He + do not coexist in comparable fractions most of the time. Moreover, the contribution of -Fpoi to the 
total free-energy remains always marginal. 
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FIG. 7: Effect of the lower boundary for <ti on the populations, the pressure, the massic entropy and the specific heat as 
a function of the density. The solid lines correspond to the case cxi > 0.8 A, the dotted lines to the case ai > 0.5 A; the 
temperature is T = 10 4 ' 7 K. 

3. Validity of the quantum correction -F qm 

As mentioned in Sec. IIII A 41 we have used the first-order term of the Wigner-Kirkwood expansion to take into 
account the quantum effects between atomic centers. This expansion becomes invalid at high density and low tem- 
perature. As a rule of thumb, the domain of reliability of the expansion is given by: log 10 Xk — log 10 p g / cm 3 ;> 2. Such 
a limitation has no consequence in an astrophysical context, as no astrophysical object with a helium composition 
exists beyond this limit. 

4- Influence of the He + -He 2+ coupling 

As mentioned in Sec. IIII D 31 the long range interaction between He 4 ", He 24 " and e~ is treated in a rather crude 
way in the present model. We have tested the influence of this approximation by submitting a few tests without this 
term. The results are illustrated in Fig. [S]for T — 10 4 5 K. We have checked other isotherms, and the conclusion is 
that the EOS and its derivatives are nearly independent of this coupling term except in a very limited temperature- 
density range. As illustrated in the next section, but also on Fig. [SJ and [SJ the reason is that pressure ionization 
occurs directly from He to He 2+ , with no regions where He 4 " and He 24 " coexist in comparable number, except at 
high temperature (T ;> 10 5 K) where temperature ionization dominates. Although we certainly cannot rule out the 
fact that this is an artefact of our model, a possible physical explanation might be the large differences between the 
ground-state energies of the different species, much larger than for hydrogen. The contribution of the ground-state 
energy of He to the total free-energy thus prevents partial ionization to occur, favoring the atomic phase. As mentioned 



14 




-2-10123 -2-10123 
log-,0 [P (g/cm 3 )] log 10 [p (g/cm 3 )] 



FIG. 8: Effect of the description of the He + -He 2+ and He + -e~ couplings on the populations (w.c = with coupling, n.c = no 
coupling), the pressure, the massic entropy and the specific heat as a function of the density. If not specified, the solid line 
corresponds to the case with coupling and the dotted line to the case without any coupling; the temperature is T = 10 4 ' 5 K. 

previously, it will be interesting to compare this prediction with experimental results and first-principle calculations, 
once they will be available, to verify whether this behaviour is a flaw of the present model or whether it reflects the 
behaviour of helium pressure ionization, an extremely interesting issue. 

B. Thermodynamical quantities 

A subset of our final EOS calculations, based on the model free-energy (I20II is presented in Tables lTllIII corresponding 
to Fig. [SI (for the pressure), Fig. ^](for the massic entropy) and Fig. ^](for the specific heat). For these calculations, 
ten internal levels have been considered, both for He and He 4 ". These ten levels are enough to represent the internal 
partition function as the highest levels are always destroyed even for the lowest density we are considering. No He 
doubly excited states have been considered in our calculation. This is a reasonable approximation because of the two 
following reasons. First of all, the high energy cost of these states (the first doubly excited state lies ~ 60 eV above the 
He ground state) disfavors their formation (in a way similar to the direct ionization of He to He 2+ without any He 4 " 
state, see following discussion). The second reason is their rapid decay by autoionization (typically in 10~ 13 — 1CP 14 s). 
It is therefore unlikely that these states survive in the midst of interacting neighbour particles. The zero of energy 
corresponds to the fully ionized plasma at zero temperature. The rising behaviour of Cy for log 10 p g / cm 3 ^> — 1 stems 
from correlations between helium atoms (configurational free-energy), since all excited levels are destroyed at this 
density, at least for the coolest temperature. The drop at larger density reflects pressure ionization, from He to He 2+ . 
We also present in Fig. ^Jthe phase diagram of helium. The lines separate the different domains where either He, He 4 " 
or He 2+ is the dominant species, i.e. represents a fraction larger than 50 %. As mentioned previously, an interesting 
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TABLE I: Equation of state for the isotherm T = 10 3 ' 8 K. For each value of the density are given the abundances of He, He + 
and He 2+ , the pressure, the massic entropy, the massic internal energy (with a zero of energy corresponding to the fully ionized 
plasma at zero temperature) and the specific heat. 
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prediction of this diagram (see also the Tables) is that for T <J 10 5 K, pressure ionization, defined as x He 2+ ;> 0.5, 
proceeds directly from He to He 2+ at p 10 g cm -3 , i.e. P > 20 Mbar. As mentioned in Sec. IIV A 41 the sharp 
transition due to pressure ionization, from ihc <; 0.5 to £ He 2+ ;> 0.5 at p ~ 10 g cm -3 (see Tables), and the persistence 
of atomic helium at high density, might reflect the large energy cost of the ground-state energies of ionized species 
(24.587 eV and 79.003 eV) to the total free-energy. Eventually, abrupt ionization occurs from He to He 2+ , unless 
temperature is high enough to unbound one of the two electrons from the helium atom. This is corroborated by the 
fact that the pressure ionization of He + (which happens if T > 10 5 K) occurs at lower densities, p > 1 g cm -3 . This 
phase diagram can be compared with the one for hydrogen [3lT |. For deuterium, the EOS is essentially the same as 
for hydrogen providing the nucleus mass is rescaled |5j. However, for helium, because of the Z = 2 nucleus and the 
induced electronic structures, the phase diagram is different, and pressure ionization occurs at larger pressures than 
for H or D. 



V. CONCLUSION 



In this paper, we have computed a free-energy model aimed at deriving the thermodynamic quantities of dense 
fluid helium, from the low-density atomic domain to the high-density fully ionized regime, covering the regime of 
partial ionization. The model is based on the so-called chemical picture for the description of the interactions between 
the different species in the fluid. The abundances of the various atomic and ionic components are obtained through 
minimization of the free-energy Despite the shortcomings inherent to the chemical approach, wc believe the present 
model to give a reasonable description of the equation of state of dense helium, including the regime of pressure 
ionization. Although the basis of the model become of doubtful validity in this latter domain, this affects only limited 
regions of the temperature-density diagram. Comparisons with available sound speed measurements and shock-wave 
experiments for atomic helium assess the validity of the model up to the megabar range, whereas at very high density 
the model recovers the fully ionized plasma model and thus Monte-Carlo simulations of the thermodynamic properties 
of the so-called one-component plasma (OCP) model. Although the present model cannot pretend giving a precise 
determination of the various atomic and ionic concentrations in the fluid, at least in the pressure ionization regimes, it 
yields a reasonably accurate determination of the phase diagram of dense, fluid helium with its various He/He + /He 2+ 
ionization contours. For T < 10 5 K, pressure ionization is found to occur directly from atomic helium He to fully 
ionized helium He 2+ , or at least to a strongly ionized state, without He + stage (cc Hc + < a few %). It would be 
interesting to test such a prediction with high-pressure dynamical experiments. Indeed, such a behaviour of the phase 
diagram bears important consequences for the thermodynamic, magnetic and transport properties of the interior 
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TABLE II: Equation of state for the isotherm T = 10 4 ' 2 K. For each value of the density are given the abundances of He, He + 
and He 2+ , the pressure, the massic entropy, the massic internal energy (with a zero of energy corresponding to the fully ionized 
plasma at zero temperature) and the specific heat. 
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TABLE III: Equation of state for the isotherm T = 10 4 ' 5 K. For each value of the density are given the abundances of He, He + 
and He 2+ , the pressure, the massic entropy, the massic internal energy (with a zero of energy corresponding to the fully ionized 
plasma at zero temperature) and the specific heat. 
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of cool and dense astrophysical objects, including giant planets. In all cases, pressure ionization is found to occur 
around p <~ 10 g cm" 3 , i.e. P <~ 20 x 10 6 bar. Detailed explorations of the sensitivity of the results upon various 
approximations entering the free-energy model show that they remain inconsequential on the first derivatives of the 
free-energy over most of the phase diagram. In some limited regions, however, characteristic of the pressure ionization 
regime, maximum variations of the entropy and the pressure can reach ~ 5% and ~ 20%, respectively, in the worst 
case. Although still modest in most cases, the uncertainties become larger for second derivatives, in particular the 
ones directly related to the different degrees of freedom and thus to the relative populations, like the specific heat. 
As mentioned above, however, only limited regions of the phase diagram are concerned by the regime where various 
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species coexist in comparable numbers. As a whole, the present model remains simple enough to allow the calculation 
of the EOS of dense helium over an extended domain of pressure and density, a necessary condition for applications 
to the computation of stellar and giant planet internal structure and high-pressure experiment diagnostics. 

Besides its astrophysical interest, the calculation of the phase diagram of dense helium is of intrinsic theoretical 
interest. Indeed, comparison bctwen these calculations and near-future high-pressure shock-wave or laser experiments 
will allow a better determination of the domains of validity of the present model and of the possible improvements. 
By such, these comparisons will yield a better understanding of the properties of matter under extreme conditions, 
and more specifically of the complex regime of matter pressure ionization and metallization. 
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